#include <dune/initial/initial.hh>
#include <dune/boundary/boundary.hh>

// initial condition
namespace Dune {
  namespace Compose {
    template<int dim, int nvar> 
    class InitialWrapper {
    public:
      typedef double RF;
      using Range = Dune::FieldVector<RF,nvar>;
      template<typename INIT>
      InitialWrapper(const INIT& init){
        this->init = init;
      }
      virtual ~InitialWrapper(){
      }
      template<typename E, typename X>
      Range initial (const E& e, const X& x) const
      {
        X xglobal = e.geometry().global(x);
        Range u(1.0);
        double xtmp[dim];
        double utmp[nvar];
        for (int i=0; i<dim; i++)
          xtmp[i]=xglobal[i]; 
        init->initial_condition(xtmp,utmp);
        for (int i=0; i<nvar; i++)
          u[i]=utmp[i]; 
//        std::cout << "call user initial condition, xglobal: " <<  
//           xglobal[0] << " " << xglobal[1] << ", utmp: " << 
//           utmp[0] << " " << utmp[1] << " " << utmp[2] << std::endl;
        return u;
      } 
    private:
      Dune::Initial::InitialBase* init;
    };

    template<int dim, int nvar> 
    class BoundaryWrapper {
    public:
      typedef double RF;      
      using Range = Dune::FieldVector<RF,nvar>;
      template<typename BOUN>
      BoundaryWrapper(BOUN& boun){
        this->boun = boun;
      }
      virtual ~BoundaryWrapper(){
      }
      template<typename E, typename X>
      Range fixed_boundary (const E& e, const X& x) const
      {
        X xglobal = e.geometry().global(x);
        Range u(1.0);
        double xtmp[dim];
        double utmp[nvar];
        for (int i=0; i<dim; i++)
          xtmp[i]=xglobal[i];
        boun->fixed_boundary(xtmp,utmp);
        for (int i=0; i<nvar; i++)
          u[i]=utmp[i];
        return u;
      }
    private:
      Dune::Boundary::Boundary* boun;
    }; 
  } // Wrapper
} // Dune
